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SUMMARY 


This report contains the results of an investigation carried out for the 
NASA Langley Research Center TCV B-737 program to study the problems in 
navigation and guidance encountered by aircraft in the initial transition period 
in changing from DME, VORTAC, and barometric instruments to the more 
precise MLS data tji)e navaids in the terminal area. The report investi- 
gates the effects of the resulting discontinuities on the estimates of position 
and velocity for both optimal (Kalman type navigation schemes) and fixed gain 
(complementary type) navigation filters, and the effects of the errors in cross 
track, track angle, and altitude on the guidance equations and control commands 
during the critical landing phase. A method is presented to remove the dis- 
continuities from the navigation loop and to reconstruct an RNAV path designed 
to land the aircraft with minimal turns and altitude changes. 

Areas for additional study are recommended to cover special problems 
encountered in the study. Smooth transition performance is shown to be pos- 
sible and feasible for aircraft utilizing MLS navaids in the terminal area. 



INTRODUCTION 


During the noncritical phases of automated aircraft flight, navaids 

such as VORTAC, TAG AN, and baro-altimeter are considered sufficiently 

accurate for navigation, guidance, and control. Range bias errors of the 

order of magnitude of 500 meters are not uncommon, heading bias errors 
o 

of 1 are often encountered, and corresponding altitude bias errors of 40 
meters are ordinarily met. With the introduction of MLS in the terminal 
area, these instrument biases are abruptly reduced to 30 meters in range, 

. 05^ in azimuth and elevation, and 5 meters in altitude. 

The effect of this sudden realization of the true error is to present 
the navigation and guidance system with a large deviation in cross track, 
track angle and altitude. These command the control system to execute 
large roll angles to wipe out the cross track and track angle errors, and to 
call for large elevator and throttle changes to attain the desired altitude and 
altitude rate. In addition, the navigation filter, in response to the large 
residuals, causes transient errors in the estimate of the position and velocity. 
Finally, the guidance gains in the terminal area tend to be somewhat higher than 
those used during the noncritical flight phase which further aggravates the 
transition problem. 

The puipose of this study is to simulate typical transitions in the 
terminal area and to investigate what software changes might be undertaken 
to reduce the severity and perhaps eliminate the navigation guidance and con- 
trol problems described above without unduly changing the existing guidance 
laws or navaids. 

Simulation programs (ALERT, BANKTURN) of the TCV B737 aircraft 

* 

are described in Reference 1. In order to simulate the transition problem it 
was necessary to add a simulation of the RNAV waypoint path design equations, 
the noncritical navaids such as VORTAC, and the baro altimeter and to duplicate 


1. Pines, S. ; Schmidt, S. F. *, and Mann, F. : Automated Landing, Rollout, 
and Turnoff Using MLS and Magnetic Cable Sensors. NASA CR-2907 . Oct. 1977. 



the guidance equations for the lateral and vertical path control. Finally, it 
was necessary to include a simulation of both the Kalman and complementary 
navigation filters in a single program in order to be able to illustrate the 
relative estimation performance of a given filter when the alternate filter was 
being used in the guidance loop. The program containing these changes is 
named FILCOMP, Only those items that are new and not contained in Ref. 1 
are described herein. The analytical work is contained in four Appendices. 
Appendix A contains an introduction into vector and matrix operations which 
are extensively used in the waypoint path equations and in the guidance for- 
mulation. Appendix B contains the derivation of the wa 3 q>oint path equations 
and guidance. Appendix C contains the derivation of the third order com- 
plementary filter in discrete form (i. e., without internal numerical integration 
loops) which is not used in the TCV B737 flight control computer. Appendix D 
contains the derivation of navigation filter equations for estimating the wind 
direction and speed. 

The main body of the report contains the results of the simulation study 
carried out using the FILCOMP program. 
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I. RNAV GUIDANCE EQUATIONS 
( 3) Coordinate Frames 

The RNAV path equations are defined and derived in Appendix B. 
The desired horizontal path consists of a series of straight lines (great 
circles) connected by circular turns of constant radius. The vertical path 
is a piecewise linear function of the ground track distance between the con- 
secutive midpoints of the circular turns. The entire path is fixed in the 
Earth and rotates rigidly with the Earth about the North Pole. The guidance 
equations are designed to null out the cross track, track angle and altitude 
errors between the aircraft position and heading with respect to this rotating 
path. In order to accomplish this, the navigation system must provide esti- 
mates of these errors from knowledge of the desired position and velocity 
and the best estimate of the aircraft actual position and velocity. 

Several coordinate systems are required to accomplish this task. 

The primary inertial coordinate system is illustrated in Fig. la where x^^ is 

directed along the Earth North Pole, x is fixed at the intersection of the 

ol 

Greenwich Meridian and the Earth equator. The x coordinate forms a 
right hand orthogonal system defined by 



] = Cx3j3 X { 


A 


X 


II 


( 1 ) 


The second coordinate frame of interest is the Earth fixed coordinate 
system. In this reference frame, x is along the Earth North Pole, x is 

-Lilt oE 

fixed along the intersection of the Greenwich Meridian and the Earth equator 

and rotates with the Earth at the earth rotation rate, w ; x forms a right 

E 2E 

hand system. The transfoimation from the E frame to the I frame is given by 


^lE = ^ < ^E 


(2a) 
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(b) Earth Fixed Inertial Coordinate System 


(GREENWICH AT t ) 


Figure 1. -Concluded. 



The transformation matrix, T (a , x ) is defined in Appendix A. The 
coordinate system is illustrated in Fig. 1(b). We have 


TE 


1 

0 

0 


0 


0 





The value of 

E 


is given in Table I. 


(2b) 


To obtain the vector velocity and acceleration of a point in Earth fixed 
coordinates in terms of inertial coordinates, we have 


C Rj, } = f ^ ■ "'e ^ ^ ^ ^ ^ 

{ Rj, ] = { Rj } + w| { Rj } (2d) 


The third coordinate frame used in the guidance equations is the local 
level coordinate frame. In this system the x coordinate points due north, 

i~±j 

the X- coordinate points due west, and the x coordinate is positive up away 

Zilj oLi 

from the Earth center through the center of gravity of the aircraft. The origin 
of the system is at the center of gravity of the aircraft. The transformation from 
local level to inertial coordinates, T , is given by 

Tj^ = T ( , x^ ) T ( «2 . Xg ) (3a) 

where 

®2 = 6 

X is the longitude of the aircraft 
6 is the latitude of the aircraft 
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The coordinate system is illustrated in Figure 2(a). The transformation, Tj^ 
is given by 


cos 6 


sin 6 


IL 


sin 6 sin a- 


-sin 6 cos a. 


cos 


sin 


-cos 6 sin 


cos 6 cos a. 


(3b) 


J 


Given the inertial coordinates of the aircraft ( , x^^ ) we 

have the unit north vector, 


r A 
cos 0 I 

\ 

^ i 

N = ( sin 6 sin \ 

l^-sin 6 cos ^ 
the unit west vector, 

f" 1 

* j i 

W = cos ^ 

I * ^ 

i^sin 0^ J 
and the unit radius vector, 

A 

R 

where 

6 



(4a) 


(4b) 


(40) 


(4d) 
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The aircraft longitude is given by 


X (t) = - Wg 


(4e) 


A fourth coordinate system is the level ground coordinate which is a rotation 
of the local level system about the north pointing axis through 180^ (see Fig. 2b). 
Thus, X is north, x is east, and x is down. We have 

ICr oCj 



T (a^ 


x^ ) T ( 6 , ) T ( 7T , ) 


( 5 ) 


To find the aircraft track angle we first obtain the relative velocity 

0 

vector with respect to the Earth, R , expressed in the inertial frame and 

REL 

then rotate this relative velocity vector into the local ground system. 


{R(,} = (T^j) 


(6a) 


where R 


REL = ^ - "'e * ^ ^ 


Finally, 




(6b) 


Ro(2) 


The aircraft altitude rate, h , with respect to the ground is given by 


h = (3) (6c) 

Similarly for the vertical acceleration in the local ground system we 

have 

{ R^ } - T^j f ^ 

^REL = ^ " ^E^ ^ 

and h = -R (3) ('^) 

(_r 
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The aircraft altitude is given by 


h = 



2 


+ X 


21 


+ 




( 8 ) 


The final coordinate system used in the guidance equations is the 

body system. The x. coordinate is along the fuselage pointing forward, 

-Lj3 

the X coordinate vertically down toward the level ground, and the x 

36 26 

coordinate is positive out along the right wing forming a right hand system. 
See Fig. 3. The transformation from body axis to inertial is given by 


(b) RNAV Guidance Parameters 

The desired path is divided into segments. Each leg is either a straight 
line (great circle) terminating at the intersection with an arc of a circle, or a 
portion of a circular turn terminating at either the second half of the same 
circular turn or at the intersection with a straight line (great circle). The 
end of each leg is defined by a fixed vector waypoint whose coordinates are 
generated in the fixed Earth coordinate system in Appendix B. Along each 
segment we are required to generate the following RNAV parameters re- 
quired for guidance: 

(1) Cross track error 

(2) Track angle error 

(3) Distance to go to the end of the leg 

(4) Time to go to the end of the leg 

( 5) Desired altitude and altitude rate 

(6) Desired airspeed 

(7) Desired flight path angle 

(8) Desired bank angle in the turn 
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A 



Figure 3(a).- Body Coordinate System Relative To Level Ground 
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The equations for these quantities are different on a straight line 
than they are along a circular turn. The equations will be derived 
separately for each segment type. We may assume that we have available 
from the Inertial Navigation IMt the position, velocity and acceleration 
vectors in the inertial coordinate system. 


(1) Great C ircle Guidance Parameters 

A 

The desired position on a great circle, WD , lies in a plane which 

A 

contains the unit radius vector to the aircraft, R_ , and the normal to the plane 

A A 

containing the great circle, WN . WN is generated for each great circle by 
means of the equations given in Appendix B. Since every line in a plane is 
perpendicular to the normal of that plane, it follows that 


WD 


1 

cos 


R 


E 


sin j3 
cos ^ 


WN 


where 


sin jS 
cos /3 


= -WN • R^ = 

= y 1 - sin^ P 


-cos ( 7 t/2 + /3) 


(10a) 


(10b) 


See Fig. 4. 


The cross track error is given in terms of the angle between R_ and 

A A A K 

WD. Since the normal, WN,is perpendicular to WD, we have 

^ A A 

CRTE = “ r_ sin"^ ( R • WN ) (11) 

£i £j 


The distance to go to the end of the s^ment is given in terms of the angle 

A 

between WD and the waypoint at the terminus of the great circle segment, 

A 

PI(J) , which is computed in Appendix B. 


DISTGO = r^ sin (| WD x PI(J) |) 


( 12 ) 
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To obtain the desired heading along the ground track at the desired 

A A 

path position vector, WD , we note that the unit vector, WD, is perpendicular 

A 

to the plane which contains the unit north vector, N , the unit west vector, 

A A 

W , and the unit normal, WN . (See Fig. 5). 

The desired heading, ^ , is the negative of the angle through which 

13 A 

we must rotate the rinit normal vector, WN., about the desired unit position 

A A 

vector, WD , to have it coincide with the unit west vector W. Thus 

A A A A 

W = cos WN - sin WD x WN 


It follows that 


j/) - tan 


WN • N 

5C“ 

WN • W 


The track angle error is defined as ^ . Since this ai^le is small 

it is more accurate to use the tai^ent of the difference of the two angles rather 
than take the difference of two arc tangents. We have 


_ tan 0 - tan 0 

TAGE = tan ^ 

1 + tan 0^ tan 0^ 


The recommended equation is given by 

R (2) (WN • W) - R^(l) (WN • N) 

TAGE = tan f ^ x x 

R^(l) (WN • W) + R-(2) (WN • N) 

vjr ta 


To obtain the time to go to the end point of the great circle segment 
we divide the distance to go by the ground speed and obtain 


TMTGO = 


DISTGO/VG 


’A 


where 


VG 


( 1 ) + 


Rg"(2) 


1 

2 


(13a) 

(13b) 

(14fl) 

(14b) 

(15a) 

(ISb) 
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The RNAV equations for desired height, HC , the desired altitude 
rate, HDC , and the desired airspeed are given by 

(16a) 
(16b) 
(16c) 


HC(t) = HEND - GRAD • DISTGO(t) 
HDC(t) = VG(t) • tan (GRAD) 
VDESm(t) = VGE - DVG * DISTGO(t) 


where, 

HEND is the desired height at the end point of the segment and is 
computed in Appendix B 

GRAD is the desired altitude gradient and is computed in Appendix B 

DVG is the desired airspeed gradient (change in airspeed per unit 
of groimd track distance) and is computed in Appendix B. 

VGE is the desired airspeed at the end of the segment aircraft is 
tracking. 

The bank angle command for the great circle is zero since we desire a 
straight line path. 


BANK = 0. 


(17a) 


In order to provide lead time to execute the bank command at the 
onset of the turn, the bank command is switched to the value it would have 
in the turn, namely. 


BANK = sign ( A0 ) tan 



(17b) 


as soon as the time to go to the end of the great circle is less than the bank 

command divided by the maximum bank rate, ^ . We switch from 

MAX 

Eq, (17a) to Eq. (17b) whenever 


TMTGO s tan 


<VGyYG/*MAX 


(17c) 
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where 


Rtji is the turn radius of the oncoming turn, and G is the constant 
of gravity. 


and 

A0 is the change in heading to be executed by the turn. A*/) and 
its sign are computed in Appendix B. 


(2) G uidance Pa rameters in the Turn 

A 

The logic for defining the unit desired vector, WD, on the turn circle 

A 

and the unit normal to the path at the turn circle, WN , is based on the con- 
cept that they both lie in the plane defined by the unit aircraft position vector, 

A A 

, and the unit center of the turn vector, CR . See Fig. 6. The vector 

^ A A 

equation for WD and WN are given by 

A ^ A A 

WD = — 3- (sin ot RE “ sin { a- j3 ) CR) 
sin jS \ f / 

and 

WN = - ■ ( - cos 0? RE + cos (c^ - ]8) CR) 

sin ^ 

where 





and 

/3 = sin ^ I RE x CR | 

Having obtained the desired unit vector position on the turn and the normal to 
the desired path, the equations for the guidance parameters for the cross 
track error, CTRE , and the track angle error, TAGE , are identical to 
those given in Eq’s. (11) and (14b). 


(18a) 


(18b) 


(18c) 




The distance to go equation to the end of the turn s^ment differs 
from the straight line DISTGK) in two respects. First the ground track 
length of the turn is the arc of a circle whose radius is and not r^; 
second, the turn is divided into two segments in order to provide a 
mechanism for introducing changes in the altitude and airspeed gradients 
which may occur at the midpoint of the turn. 

To solve for the distance to go on a turn we compute the angle between 

A 

the unit normal at WD and the unit normal at the end of the circular segment, 

A 

WNE. The normal at the end of the segment is computed in Appendix B. We 
have 


^ A A 

e = sin"^ ( I WN X WNE | ) 

where WNE is the normal at the end of the circular segment 

Since the sin 6 = sin ( - 6 ) we must test to see if 6 < /2. To accom- 

plish this we compute 

A A 

COS 0 = WN • WNE 

If cos 6 is negative, we set 

0 = 77 - Sin” ( I WN X WNE | ) 

The final result yields 

DISTGO = • 0 

The equations for the time to the end of the segment, TMTGO , HC(t), 

HDC(t) and VDESIR(t) are all computed in a manner identical to Eq’s, (15a), 
(15b), (16a), (16b), and (16c). 


(19a) 


(19b) 


(19c) 


(20) 
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The bank angle command logic aloi^ the turn is given by 


BANK = sign (A0 ) tan ^ (VG^yo/R^ ) (2X) 

The lead command for anticipated return to the straight line is exer- 
cised only at the end of the second half of the turn. The bank angle command 
is continuous across the middle of the turn. Toward the end of the second half 
turn we test to see if the time to go is less than the time to return the bank angle 
to zero at the maximiun bank rate. Whenever 

TMTGO s tan'^ 

we switch the bank command from Eq. (21) to BANK - 0 . 

(3) Switching Logic 

The segments are arranged in groups of three. A great circle followed 
by a half turn, followed by another half turn. Only the final leg terminating at 
a touchdown waypoint is treated separately. Each leg is numbered in sequence 
with the first leg initialized as segment number 1. This is always a straight 
line segment. The guidance information common to a given segment, such as 
the unit end waypoint, the desired height, airspeed and gradient, etc. are 
stored in an array for that segment. These data are computed once for the 
entire flight following the equations outlined in Appendix B. 

If the segment number, modulo 3, is 1 we are on a great circle straight 
line segment and enter the data block for active guidance computations for a 
great circle arc. This arc is terminated whenever 

TMTGO < .075 seconds (23) 
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A logical trigger defining a great circle arc as being TRUE is set 
to FALSE and remains FALSE imtil the segment number, modulo 3, is 
equal to 0, The segment number is advanced by unity at the end of each 
segment and the block storage for the segment parameters is brought into the 
active guidance loop. 

On the final leg a test is made on the last segment number and if the 
leg segment is equal to the last integer no further call to the block storage 
is made. 
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II. 


NAVIGATION AND FILTERING FOR WAYPOINT GUIDANCE 


(a) T he Estimate of _the Inertial State 

The aircraft is equipped with an inertial navigation system (INS) 
consisting of body mounted accelerometers for measuring the specific 
forces acting on the aircraft, plus a stabilized platform containing rate 
gyros for integration of the aircraft attitude and attitude rates. In the 
absence of instrument errors the inertial position and velocity may be 
obtained by integrating the inertial equations of motion 


dt r^ 

where 

^ 2 
u = G 

^ E 

rj =KI 

A A A A A A 

Tib = T( 0 ^ 1 , X^) T(fl^,X 2 ) T(ff,x^) T( 0 ,X 3 ) T( 0 ,X 2 ) T(( 0 ,Xi) 


(24a) 


(24b) 

(24c) 

(24d) 


{f } is the specific force vector in the body system as computed 
SFB 

by the IMU from the output of the body moimted accelerometers, [ 



^ *SFB ^ 


cc 


(24e) 



(24f) 
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The Euler angles <p , 6 , and 0 are computed by the IMU system from the 
output of the rate gyros mounted on the stable platform. 

Since the rate gyros are subject to drift, the accelerometers contain 
biases and all the measurements are noisy, it is required to design an esti- 
mator to filter out the noise and obtain the best estimate of the aircraft state 
and the instrument biases. In addition, external navaids are required, such as 
VORTAC range and bearing as well as an aircraft barometer during the non- 
critical flight phase. For the critical terminal area, MLS range, azimuth, 
and elevation are used. These navaids are themselves subject to instrument 

biases. 

(b) Filter State Estimators 

Two estimators are described here. First, an optimal Kalman filter 
and second, a fixed gain complementary filter. Both are used in the study in 
a manner that permits the aircraft guidance equations to utilize the output of 
one filter with the non-controlling filter operating only in its estimating mode. 
In this manner both filters provide state estimates for comparison with the 
true state, and both filters can be compared as guidance state estimators as 
well. 

The two filters are designed to operate in a local flat Earth coordinate 
frame with the origin set at the navaid site. Thus when the VORTAC mode 
is in operation, the coordinate system for the error state estimate has its 
origin at the VORTAC station site 

A 

X, is due north 
Iv 

A 

x^ is due east 
2v 

A 

X is down toward the Earth center 
3v 
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The coordinates of the VORTAC station in the inertial system are 


given in terms of the station longitude, latitude , and altitude, h^ , 

sin 6^ 

C Rj^(t) } = f sin cos 6^, ^ ) 


IE 


E 


cos Xy cos 6y 


( 25 ) 


To obtain the local flat E&rth coordinates of the aircraft with respect to the 
VORTAC station, we have 


^ ^FV ^ ^FI ^ ^ ^ 


(26a) 


where 


T 

FI 

A A A 

= T( - ff, x^) T ( - ttg, x^) T ( -ftj, x^) 

(26b) 

r-l 

= wjj(t-y +Xy 

(26 c) 

“2 

II 

o» 

<1 

(26d) 


( See Fig. 7 ) 


In the case of the MLS landii^ site, the local flat Earth coordinate 

A 

system is oriented so that is aligned with the runway and the origin is 
located at the center line of the runway as shown in Fig. 8 . 

The inertial coordinates of the runway origin are given by 


{ } 


TE 




cos 

cos 



(T 


E 


h ) 

r' 


(27) 
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Figure 7. - VORTAC Flat Earth Coordinate System 
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Figure 8. -Runway Flat Earth Coordinate System 



To obtain the local flat Earth coordinates of the aircraft in the runway system, 
we have 


f = ^FI 


(28a) 


where 


Tpj = T( - 0^, Xg) T ( -ff, x^) T ( -Og, x^) T ( x^) 


(28b) 


\j) - runway azimuth measured from true north 

R 


“2 = «R 


(28c ) 


The Kalman error state is a 15 element vector consisting of 


R 


error in the estimate of the aircraft position vector in 
the flat Earth coordinate system 


R 


P . q , r 


w 

X 


> 


w 

y 


error in the aircraft estimate of the aircraft velocity 
vector in the local flat Earth system 

error in the estimate of gyro drift rates 

error in the estimates of horizontal components of the 
local wind 


b 

b 


1 

4 


» 



errors in the estimate of the navaid instrument biases 

error in the estimate of vertical acceleration due to 
instrument bias . 


When the VORTAC navaids are operative, we have 
b^ = bias in DME (range) 

bg = bias in bearing 


bias in the baro reading of pressure altitude 



Wheh the MLS navaids are operative, we have 
= bias in azimuth 

b = bias in elevation 

b = bias in range (DME) 

O 

Immediately prior to landing, at an altitude of 50 meters, the elevation reading 
is replaced by a vertical radar and b is the estimate in the bias of radar 
altitude. 

The filter equations for the Kalman optimal estimator for a flat Earth 
coordinate system are described in Ref. 1 and are not repeated here. 

The error state for the fixed gain complementary filter consists of 
the following; 

R Three components in the error of the position vector in the 
flat Earth coordinate system 

R Three components in the error of the velocity vector in the flat 
Earth coordinate system 

• • 

R Three components in the error of the estimate of acceleration 
in the flat Earth coordinate system 

w^,Wy errors in the estimates of the horizontal components of the winds 

Thus the complementary filter is a nine state filter for the aircraft position, 
velocity, and acceleration. This is amply described in Ref. 1 and is not 
repeated here. The derivation of the complementary filter in discrete 
explicit form is derived in Appendix C since the derivation was not included 
in Ref. 1. 
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The navigation equations for the two filters are different and are 
described herein, 

(1) Kalman Navigation Equations 

When using the optimal Kalman filter, the navigation system estimates 
the expected bias errors in the gyro drifts. These are converted into correc- 
tions in the Euler angles and added to the output of the IMU estimate of the 
Euler angles as described in Ref. 1. 

A 

® = ®IMU 

The Kalman filter also estimates the bias error in the vertical accelerometer. 
This is treated as a correction to the estimate of the vertical specific force 
in the body system. 



^SFB ^SFB - ^4 

(30a) 

and 

^SFB ~ ^SFB 

(30b) 


A 

f (2) = f (2) 

SFB ' ^ SFB ' ' 

(30c) 


The Kalman navigation equations integrate Eq. (24a) in the inertial reference 
frame using an Euler integration foimula. 

A A A ^ 

• • A.t * * * ■ 

{ R2 (t) } = C Rj (t - At) } + -^ { Rj (t) + Rj (t - At) ] (31a) 
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where 


A, 

A A 

c ii; (t) ) = - — [ IL (t - At) } + T [ f 3 

r (t-At) ^ 

and the angles used in computing the transformation T are 

IB 

01 ^ = tan ^ - At) / Kj 3 ( t - At)^ 

6^ = sin“^ (t - At) / (t - At)) 

A 

j/3 = 0 (t) 

A 

6 = d (t) 

A 

<P = <p (t) 


(31b) 


(31c) 


The Kalman filter velocity correction is "eased in" in an exponential manner so 
as not to present the control command with abrupt discontinuities as described 
in Ref. 1. 


[ Rj (t) } = C Rj (t) 3 


+ e *k [ R (t^) 3 ) 


(32) 


where the angles used in computing the transformation Tjp , are given by Eq’s. 

(26b) - (26d) or Eq’s. (28b), (28c) depending on the station coordinates. 

The time, t^ , is the last Kalman update time prior to the present time, 
t . The Kalman filter equations are processed only once every 10 normal time 
steps as described in Ref. 1. 

The integration of the velocity in the inertial reference frame to obtain 
the position vector is accomplished in a similar fashion. 

[ Rj (t) 3 = f Rj (t) 3 + ^ C Rj (t) + Rj (t - At) 3 (33) 
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The Kalman correction to the position is given by 


{Rj(t)} = {Rj(t)3 +e •(Tjj,{R(y)} 


The best estimate of the inertial state is presented to the RNAV Guidance 
and Display routines. A block diagram of the Kalman filter navigation logic is 
shown in Fig. 9. 
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Figure 9.— Navigation Equations for Kalman Filter 













(2) Complementary Filter Navigation Equations 


The complementary filter navigation equations are described in Appendix C, 

A prediction of the aircraft position in inertial coordinates is made using the 
assumption that the acceleration is constant over the start time interval, At. 

The Euler angles, (p , 0 and 0 are accepted without correction directly from 
the IMU. We have 

A A 

{ E^(t) } = ' { Rj(t-At) ) + Tjg c fgpgW] + T C R (t-At) } (35a) 

rj(t-At) 

the angles used in computing when the transformation T are 

IB 

= tan C- ^3 ~ 

«2 = 

(p , 0 and 0 

and angles for computing the transformation T are either Eq’s. (26b) - (26 d), 

IF 

or Eq’s. (28b), (28c) depending on the station coordinates. The vector 
[ R ( t - At) } is the previous estimate of the acceleration random error 
variable from the complementaiy filter. 

To obtain the uncorrected velocity and position in inertial coordinates, 
we have 

A A A 

[ Rj(t) } = [ Rj (t-At) 3 + At { 'Rj(t) 3 (36) 

A 2 ^ 

{ Rj(t) 3 = [ Rj (t - At)3 + At C Rj(t - At)3 + [ i^(t) 3 (37) 

The complementaiy filter prediction of the pseudo observation used in the 
complementary filter is the position of the aircraft relative to the navaid site in 
the flat Earth coordinate system given by either Eq. (26 a) or Eq, (28a). 
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We have 

{ Rpg ) = Tpj c ^ (t) - Rjg (t) 3 (38) 

when { Rjg (t) } is the inertial position vector of the navaid site. 

The computation of the pseudo observation of the relative position in 
terms of VORTAC range, r^ , VORTAC bearing , and bare altitude, 

, is given by 

Rj;^(l) = cos ( 0^) (39a) 

R^„(2) = r sin ( 0 ) 

FS V ^v 

- 2 2 

= "Z +(r -z )/2r^+h 

Eq. (39c) yields the height above the plane, tangent to the Earth at navaid site 
as a function of the altitude above the spherical Earth and the range to the 
navaid site. 

The computation of the pseudo observation in terms of the MLS i^nge, 

^MLS* MLS azimuth, elevation, is given in Ref. 1. 

The complementary filter estimates of the errors in position, velocity, 
and accelerations are 

f - ^FS 5 

M «S» A 

= ‘^FS<* - ^ S f ^FS<*> - «FS ’ 


(40a) 

(40b) 

(40c) 


(39b) 

(39c) 
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The best estimate of the relative state is given by 


+ { ^g(t)} 




( 41 ) 


To obtain the inertial state required for ^idance and display routines, 
we have 


A 


CRl(t)} 


(42a) 

A 

f^(t)} 

'' A A 

“ ’’if f '^E 

(42b) 


A ^ 

(42c) 


A block diagram of the complementary filter navigation logic is presented 
in Fig. 10. 
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Figure 10.- Navigation Equations for Complementary Filter 
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NAVIGATION AND GUIDANCE AT THE TRANSITION POINT 


(a) Filter Update at Transition 

Both the optimal Kalman filter and the fixed gain complementary filter 
require a settling time period when presented with an abrupt change in navaid 
observation type. As previously discussed, during the transition from VORTAC 
to the more accurate MLS data, the prediction of the MLS observations, based 
on the estimates of the aircraft position derived from VORTAC measurements, 
result in a large residual. In order to avoid undesirable effects of this pertur- 
bation , it is most convenient to reinitialize the aircraft estimate of its position 
from the first valid MLS observation set. The equations necessary to accomplish 
the change are presented for both filters. 

(1) Kalman Filter Update at Transition 

Let the coordinates of the range and azimuth antenna, in runway flat 

Earth coordinates be, f Raz } » coordinates of the elevation antenna 

relative to the azimuth antenna be { R }. The position of the aircraft in the 

EL 

runway flat Earth coordinate system computed from the MLS range, rjyj , the 
MLS azimuth , , and the MLS elevation , , is given by 


"M 


{R^S? = 


AZ^ 


X 


iM 

'2M 

3M 


(43a) 


where 




^2M = 


^3M = 0 ^ < =‘2M - *2ElO 


(43b) 


^ = ^lEL < ^M> 


2 2 


2 . . 2 


^2 * “ ^M ^ ^ ^2M " ^ ^^2EL ^2M ” ^2EL ^lEL ^ ^ ^M^ 
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Prior to processing the first MLS observation, the Kalman best 
estimate of the aircraft inertial position is replaced by 




(44a) 


The angles used to generate the transformation T are the negative 

IF 

of the angles listed in Eq. (28 c) where the transformation T is defined. 

FI 

In addition, the square root of the state covariance matrix, W(t), is 
reinitialized to diagonal form in order to eliminate undesirable cross corre- 
lations that were accumulated during the VORTAC data processing. We have 
for the diagonal elements of 


15 i 

W.. = ( S W.. ) 

• 1 1 ] 

3=1 


W.. ( j db i ) = 0 
11 


and 


i = 1, 9 
i = 13 


( 443 ) 


The covariance elements of the expected error in the horizontal wind 
estimate, W. . - . , W , W and W are retained unaltered. 

The elements of covariance matrix for the VORTAC instrument biases 
are eliminated and replaced by the diagonal elements of the MLS instrument 
biases 


*10,10 = (30 meters) 

* 11,11 = 

* 12,12 = (.05°) 


(45) 
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( 2 ) Complementary Filte r Upda te at Transition 

In a manner completely similar to the Kalman filter transition 
Eq’s. (43a), (43b) and (44a), the estimate of the inertial aircraft position is 
given by 


In this manner the filter equations are no longer presented with abrupt 
perturbations and the processii^ of data is smooth. The aircraft guidance 
equations are now presented with an accurate estimate of the large cross track 
error, track angle error , and altitude error. 

(b) Path Reconstruction at Transition 

If the desired path, or the guidance equations remain unaltered at tran- 
sition, the bank angle command will call for a larger bank correction to recapture 
the desired path, and the stabilizer and thrust commands will call for large ver- 
tical path changes to achieve the desired altitude. In order to avoid these 
unnecessary maneuvers the method chosen in the study to solve the transition 
problem is to reconstruct the remainder of the desired path to touch down by 
eliminating the cross track and altitude errors. This is accomplished by simply 
accepting the updated valid MLS position as the new first desired way point and 
retaining the remainder of the data for the initial N way points that have not yet 
been encountered. By executing the path construction equations, outlined in 
Appendix B, with this new reduced data set, a path will be constructed which will 
alter only the following segments of the original path: 

( 1) The great circle segment from the present MLS position to the incoming 
tangent way point to the new next turn circle. 

( 2) The turn segment to the middle of the new turn. 

( 3) The outgoing segment from the middle of the new turn to the unchanged 
next great circle. 

( 4) The vertical path for altitude and airspeed gradients along the three 
altered horizontal segments. 


(44a) 
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The consequence of these changes is to eliminate the initial cross 
track error and to eliminate the error in the desired height at transition. 

No attempt to zero out the track angle error was undertaken in this study. 

In addition, in order to reduce the computing work load for an aircraft 
computer, a trigger was introduced to avoid the path change in the event the 
cross track error is less than a prescribed amount (say 100 meters). 

A 

Let R be the unit aircraft position vector in the fixed Earth 
£ 

coordinate frame. Then the longitude and latitude of the new first way point 
is 

(47a) 
(47b) 

Let the magnitude of the inertial position vector, , be r^ , then the 
first desired altitude is given by 


X. (1) = tan‘^ (- E^(2)/R^(3) ) 
6 (1) = sln"^ (Re<3)) 


h(l) = rj - 

and the new first desired airspeed is given by 


(47c) 


Vd(1) 


V . ^ (INS) 

airspeed 


(47d) 


Let P be integer corresponding to the last middle of the turn, WC(P), 
processed at the time of the MLS transition update. The new number of initial 
way points to be used in constructing the path is given by 


M = N - P 


(48) 
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We have for the remainder of the input data 


X(I) =X(P+I) 

6(1) = 6(P +1) 

h(I) = h(P + 1) 

Vj^(I) = v^(P + I) 

(I) = 2,M 

The radius of the turn is given by 
R^(I) = (P - 1 + 1) 

(I) =1, M-2 


Fig. 11 illustrates a typical reconstructed horizontal path. 
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Figure 11. —Reconstructed Path at Transition 
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IV. 


RESULTS OF THE SIMULATION STUDY AT TRANSITION 



This section contains the results of the computer runs of the FILCOMP 
program in the format of plots. Each run consists of three sets of plots. The 
first contains the input data of the run including the way point path generation 
information, the magnitude and direction of winds, the filter t 5 rpe used in the 
navigation, guidance, and control equations, the MLS boundary limits defining 
the transition point, and a plot of the aircraft true ground track illustrating the 
reconstruction of the path on transition. The second figure contains seven plots 
of aircraft performance in the following order: 

(1) Glide path deviation in meters for both the true deviation and the 
estimated deviation for the filter in use in the aircraft control. 

(2) Aircraft pitch angle in degrees for the true pitch, the measured 
pitch output of the IMU (used by the complementary filters) and the measured 
pitch corrected for the estimated gyro drift bias (used by the Kalman filter). 

(3) Aircraft altitude rate, measured in meters per second, for both 
the true rate of climb and the estimated rate obtained by the filter supplying the 
control equations. 

(4) Error in the estimate of the x (forward) position in runway 
coordinates, measured in meters, for both the Kalman and complementary filter. 

A 

(5) Error in the estimate of the x (lateral) position in runway 

2R 

coordinates, measured in meters, for both the Kalman and complementary filter. 

A 

(6) Error in the estimate of the x (vertical) position in runway 
coordinates, measured in meters, for both the Kalman and complementary filter. 

(7) Error in the estimate of the forward velocity ( x ) in runway 

IR 

coordinates, measured in meters per second, for both the Kalman and com- 
plementary filter. 

The third figure contains eight plots of an aircraft performance in the 
following order: 

(1) Cross track error, measured in meters, for both the true CRTE 
and the estimated CRTE for the filter used in the control. 
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(2) Track angle error, converted from degrees to time rate of 
change of cross track error by multiplying by the ground speed. Both the 
true track angle error and the estimated track angle error used by the 
guidance equations are shown, measured in meters per second. 

(3) Aircraft roll angle, measured in degrees for the true roll ar^le, 
the measured roll angle supplied by the IMU (used in the complementary filter), 
and the measured roll ai^le corrected for the gyro drift bias (used in the Kalman 
filter). 

(4) Error in the estimate of the north component of the wind, in meters 
per second, for both Kalman and complementary filter. 

(5) Error in the estimate of the west component of the wind, in meters 
per second for both the Kalman and complementary filter. 

( 6) Difference between the true desired airspeed and the true airspeed. 

A second plot also shows the difference between the true ground speed and the 
true airspeed. The curves are mirror images of one another if there are no 
winds, and differ in the presence of winds. 

(7) Error in the estimate of the lateral velocity, x , in meters per 

2R. 

second for both the Kalman and complementary filter. 

(8) Error in the estimate of the vertical velocity, x , in meters per 

oR 

second for both the Kalman and complementary filter. 

(b) Discussion of the Results 

The first case. Fig’s. (12a), (12b) and (12c), illustrates the performance 
of the aircraft, following transition, where the VORTAC and baro bias errors 
have been removed from the estimate of the position through reinitialization, but 
no path reconstruction has been attempted. Plots 4, 5, and 6 of Fig. (12b) show 
the effect of the reinitialization at 33 seconds on the errors in the estimate of 
position. Plot 1 of Fig. (12c) shows the laige CRTE presented to the guidance 
equations with no path reconstruction. Plot 3 of Fig. (12c) shows the lai^e roll 
angle (20^) at 40 seconds and the overshoot in the opposite direction of 10° at 
55 seconds as the aircraft attempts to recapture the original path. Plots 2 and 3 
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of Fig. (12b) illustrate the excessive pitch change and "roller coaster" effect 
in the rate of climb that result when the desired vertical path has not been 
reconstructed. Fig. (12a) illustrates the large shift in ground track as the 
aircraft attempts to recapture the original path. 

Case 2, Fig’s. (13a), (13b) and (13c), illustrates the perfect transition 
that would be possible if the actual vehicle position and velocity information 
were available to the guidance equations. Since the VORTAC bias errors do 
not appear in the estimate of the aircraft state, the cross track error, track 
angle error, and altitude errors are negligible both prior to and following tran- 
sition and a smooth trajectory results. 

Case 3, Fig’s. (14a), (14b) and (14c), illustrates the smooth performance 
at transition, with path reconstruction when the complementary filter supplies 
the aircraft position and velocity to the guidance equations. Examination of 
Plot 3, Fig. (14c) shows that no roll is called for during transition indicating 
that both the cross track error and track angle errors are small at transition 
(33 seconds). Similarly, Plot 2 and 3, Fig. (14b), show very small activity 
in the vertical channel as compared with the same plots for Case 1. A vast 
improvement in aircraft performance with path reconstruction is plainly in- 
dicated. Furthermore, a comparison of Fig. (13a) and (14a) shows that the 
change in ground path is minimal. 

Case 4, Fig's. (15a), (15b) and (15c), illustrates similar acceptable tran- 
sition performance when the Kalman filter is used to supply the guidance 
equations in conjunction with path reconstruction. The major difference between 
the Kalman and complementaiy filter effect on transition performance is illus- 
trated in Plots 2, and 3 of Fig’s. (14c) and (15c). The Kalman filter causes a 
slight aircraft roll (Plot 2) transition due to a non zero track angle error (Plot 3) 
caused by the error in the Kalman estimate of the lateral velocity. Examination 
of Plot 7, Fig. (15c), indicates that the Kalman estimate of lateral velocity con- 
tinues to grow after 25 secs. The complementaiy filter, lateral velocity error, 
decreases from 3. 5 m/sec at 25 seconds to 2 m/sec at transition (33 sec). 



while the Kalman estimate continues to grow to 4. 5 m/sec at transition. Fol- 
lowing transition, the Kalman estimates of velocity are much better than the 
complementary filter. Some additional study is recommended to remove this 
anomaly since it is expected that an optimal filter should be expected to out 
perform a fixed gain filter. In any case the performaice of both filters with 
path reconstruction does eliminate the objectional features of Case 1 with no 
path reconstruction. 

Case 5 illustrates the effect of winds on the transition performance 
with path reconstruction. A 10 knot wind at 165^ heading is simulated in this 
run. Comparing this case with Case 3, in which no winds were simulated, 
shows that the performance at transition, using the complementary filter is 
unaffected by winds. Even the error in the complementary filter estimate of 
the winds shown in Plots 4 and 5 of Fig. (16c) are similar to the error in the 
complementary filter estimate of the winds in Plots 4 and 5 of Fig. (14c) for 
no winds. 

Case 6 illustrates the effect of a 10 knot wind at -15° heading on the 
transition performance when the Kalman filter information is used in the 
guidance equations. Case 6 is to be compared to Case 4 in which no winds were 
acting and the Kalman filter was used in the guidance equations. The performance 
at transition is seen to be unaffected by winds for both Kalman and complementary 
filter estimates and both filters produce smooth transition using the path recon- 
struction method. 

Case 7 illustrates the effect of changing the transition point from a 60° to 
40° azimuth boimdary for the MLS elevation signals. Case 7 is to 
be compared with Case 5. Both cases differ only in the time and location on the 
original trajectory that the transition is initiated due to the 40° limit on the MLS 
elevation anteima. Comparing Fig. (18b) with Fig. (16b) and Fig. (18c) with 
Fig. (16c) it is apparent that except for the change in the time of transition, the 
performance is identical. 

o 

Case 8 illustrates the effect of a 40 azimuth boundary for the MLS elevation 
on the performance of Kalman filter. Once again the error in lateral velocity at the 
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instant of transition is seen to produce a slight track angle error following 
transition. Comparing Plot 3, Fig. (19c), with Plot 3, Fig. (17c) shows 
that the slight perturbation in track angle error following transition is moved 
in time and that smooth transition is possible for both filters and is imaffected 
by the sweep limits of the MLS elevation antenna. 

Case 9 illustrates what happens when the cross track error at transition 
is not large enough to trigger the path reconstruction. By changing the heading 
of the incoming leg it is possible to approach the transition point in a direction 
at light angles to the VORTAC range bias error. In this case the cross track 
error at the first valid MLS data point is too small to trigger the path recon- 
struction logic. The vertical error due to the baro bias and along-track error 
due to range bias now remains as a desired altitude error in the guidance 
equations and the "roller coaster" search for the proper altitude and altitude 
rate reoccurs. Plot 3 of Fig. (20b) illustrates the overshoot in altitude and 
rate of climb at transition. In order to avoid this effect it is deemed necessary 
to remove the CRTE limit at transition and cause the path redesign logic to 
occur at transition independent of the magnitude of the cross track error. 

Case 10 illustrates a 4 way point trajectory with two turns. Transition 
occurs after the first turn. The significant feature in this run is the persistence 
of the track angle error immediately following the transition point (125 secs. ). 
This is illustrated by the roll angle excursion from 125 seconds to 140 seconds 
in Plot 3, Fig. (21c). It appears that an alternate method in which both the 
cross track error and the track angle error at transition are nulled out should 
be investigated. 

Case 11 is an illustration of a 6 way point with 4 turns. The entire 
trajectory takes place within the MLS boundary so that no transition is en- 
countered. It is included here only for the purpose of illustrating the versatility 
of the way point path construction. 
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(c) Conclusions and Recommendations 

The main result of this study is that a smooth transition to MLS navaids 
in the terminal area is feasible and possible through the use of path reconstruction 
and filter reinitialization. The method, investigated in the study, of nulling the 
cross track error and altitude error with the first valid MLS aircraft position 
determination is simple and requires very little change in the existing way point 
trajectory construction software presently existent in the Langley Field TCV B-737 
aircraft navigation computer. 

The followii^ recommendations are made: 

(1) Investigate alternate methods of path reconstruction to null out 
track angle error as well as cross track and altitude error at transition. 

(2) Study the behavior of the Kalman filter at transition to obtain 
better velocity behavior in the VORTAC navaid area. 

(3) Perform path reconstruction when either the cross-track or altitude 
error exceed specified limits or j*ust always update to eliminate the code required 
for limit checking. 

(4) Investigate the changes that must be made in the existing aircraft 
software to accommodate the path reconstruction method outlined in this study and 
flight test the procedure for pilot comments and experience in transition in the 
terminal area, 

(5) Develop a method for path reconstruction for transition occurring 
during a turn. 


51 



TABLE I 


CONSTANTS USED IN WAYPOINT TRAJECTORY CONSTRUCTION 


w„ (EARTH ROTATION RATE) 
E 

G (GRAVITY CONSTANT) 

r (RADIUS OF THE EARTH) 
E 

1 knot = . 5144434 m/ sec 

1 radian = 57.2957795° 


. 7292131 10 ^ rad/sec 
/ 2 

9. 8066 m/sec 

6378.156 km 



TABLE n 


INPUT DATA FOR WAYPOINT TRAJECTORY CONSTRUCTION 
CASES 1 THROUGH 8 


I 

X (I) 

6 (I) 

h(I) 

Vj,(I) 


DEG 

DEG 

m 

m/sec 

1 

-77. 1315516 

40.2689558 

1084. 11 

74. 59 

2 

-77. 0847124 

40. 1709553 

435.16 

68.31 

3 

-77. 0232052 

40.2523725 

0 

68.31 


R^ = 2286 meters 
*^RUNWAY “ 
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X(to) " -77.1253474' 

6 (to) * 40.2615871° 
h (to) » 981. 79 m 
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ELBOUND * 50O 
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Aircraft Ground Track. 


Figure 12(a), -Case 1 MLS Transition Without Path Reconstruction 
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Figure 13(a). - Case 2 MUS Transition With Path Reconstruction 
and Perfect Navigation Path 
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Figure 13 <b). - Case 2 Continued 
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Figure 14(a). - Case 3 MLS Transition With Path Reconstruction 
With Kalman Filter Navigation 
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Figure 15(a). - Case 4 MLS Transition With Path Reconstruction 
With Kalman Filter Navigation 












Figure 15(c)* " Case 4 
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Figure 18(c). - Case 7 
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Aircraft Ground Track. 

Figure 19(a). -» Case 8 MLS Transition With Path Reconstruction 
dose to Turn 
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Si04‘U3. IQX 


O True 





Figure 19(b). « Case 8 Continued 









TABLE m 


INPUT DATA FOR WAYPOINT TRAJECTORY CONSTRUCTION 

CASE 9 


N = 3 


I 

X (I) 

6(1) 

h(I) 



DEG 

DEG 

m 

m/sec 

1 

-77. 0742769 

40. 2101441 

873.25 

74.59 

2 

-77. 1152352 

40. 1298037 

331. 56 

64.31 

3 

-77. 0232052 

40.2523725 

0 

64.31 


= 2286 meters 
’^RUNWAY ~ 
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X (to) * -77.07946582° 

6 (to) = 40.29405513° 

h (to) * 740.3 m 

WIND = 10 knot/16 5° 

AZBOUND » 60° 

ELBOUND = 60° 
COMPLEMENTARY FILTER 
PATH RECONSTRUCTION 



* 


Aircraft Ground Track 

Figure 20(a), - Case 9 MLS Transition When Cross Track Error 
Does Not Trigger Path Reconstruction 
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YfOS ER>n 






Figure 20(c), « Case 9 











TABLE rv 


INPUT DATA FOR WAYPOINT TRAJECTORY CONSTRUCTION 

CASE 10 


N = 4 


I 

X a) 

6 (I) 

h (I) 

Vjjd) 


DEG 

DEG 

m 

m/sec 

1 

-77.1276539 

40. 1986043 

643. 30 

74.59 

2 

-77. 1384963 

40.25056842 

643.87 

69.45 

3 

-77.0530063 

40.21295326 

241. 46 

64.31 

4 

-77.02320524 

40.2523725 

0 

64.31 


R^ (1) = 2286 meters 
R^ (2) = 2286 meters 

^RUNWAY " 
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X (to) * -77.16424184 
6 (to) =* 40.21641321^ 
h (to) = 639.62 m 
WIND * 0. 

AZBOUND * 60° 
ELBOUND = 40° 
KALMAN FILTER 
PATH RECONSTRUCTION 


♦ 


Aircraft Ground Track 

Figure 21(a). - Case 10 MLS Transition With Path Reconstruction 
and Multiple Turn Capability 



n)T CK.t1/5 





Q Selected Est 


O Kalm Est 
Q Comp Est 


O Kalm Est 
Q Comp Est 


O Kalm Est 
Q Comp Est 


O Kalm Est 
Q Comp Est 


Figure 21(b). - Case 10 Continued 














TABLE V 


INPUT DATA FOR WAYPOINT TRAJECTORY CONSTRUCTION 

CASE 11 


N = 6 


I 

X (I) 

6 (I) 

h(I) 

VdW 


DEG 

DEG 

m 

m/sec 

1 

-77. 04768972 

40.2599459 

1720.29 

77.16 

2 

-77.09841711 

40.192883238 

989.38 

77.16 

3 

-77. 05699524 

40. 16107542 

650. 11 

77. 16 

4 

-77.02892824 

40.19823274 

378.26 

72.02 

5 

-77. 04632437 

40. 22179639 

202.60 

64. 31 

6 

-77.02320521 

40. 25237254 

0 

64.31 


IVj, (1) = 1828.8 meters 
R^ (2) = 2286 meters 
R^ (3) = 1524 meters 

^RUNWAY “ 
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X (to) *-^7.08175208^ 

6 (to) * 40.2149138° 
h (to) 1216, 24 m 
WIND * 0. 

AZBOUND - 60° 
ELBOUND * 60° 
KALMAN FILTER 
PATH RECONSTRUCTION 



Figure 22(a), - Case 11 MLS Transition With Path Reconstruction 
and Multiple Turn Capability 
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TABLE VI 


VORTAC AND BARO SIMULATION DATA 
VORTAC SIMULATION DATA 


\ = 

-77.164894 

6 = 

V 

40. 403016* 

h ^ 

V 

45. 72 m 


VORTAC BIAS 

Range bias = 365.76 m 
Bearing bias = . 6° 

VORTAC RANDOM NOISE 
a = 48. 77 m 

V 



BAROMETER SIMULATION DATA 

BARO BIAS = 12.19 m 

a, = 1. 524 m 
bare 
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APPENDIX A 


Vectors and Matrices 

This appendix derives the equation for rotating a given vector 
about a fixed unit vector, through a given angle in Cartesian three space. 
The resulting transformation is used to describe the rotation of one Car- 
tesian coordinate system into another. Rules for vector addition, dot 
products , and cross products are reviewed. Finally, a vector matrix equa- 
tion is derived which describes the instantaneous rotation rate vector in 
the transformed coordinate system in terms of the rotation rates in the 
component systems comprising the transformation. 


A. Elementary Operations on a Vector 


A vector is an array of quantities. In this appendix we restrict our- 
selves to a Cartesian three space. Given a point in the coordinate system 

( X , X , X ), the vector, R , is most easily visualized as the directed line 
X ^ o 

segment from the origin (0. , 0. , 0. ) to the point. Thus the vector R , 
arranged as a column, is given by 


R 



(Al) 


The sum of two vectors, A + B, is a vector C , 
are given by the sum of the corresponding components of 


whose components 
A and B. 


C 




(A2) 


The dot product of two vectors, A • B , is a scalar functional whose 
value is given by the sum of the products of the corresponding elements. The 
operation is commutative 


A • B = B - A = b^ 


(A3) 
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Thus, A * A is the square of the Cartesian length of the vector; 


* A 2 2 2 2 

A • A = a + a + a = a 

X A o 


(A4a) 


and we have for the magnitude of A , 


a = A = (A • Ay 


(A4b) 


The unit vector , A , is given by 


A = — A 
a 


(A5) 


If we consider the dot product of the sum of two vectors into itself, 


we have 


(A+B) - (A+B) = a^+b^+2A*B 

From the cosine law in trigonometry, we have, if C = A + B 
2 2 2 

c = a+b+2ab cos 0 
Thus, if 0 is the angle between A and B, we have 


(A6) 


(A7a) 


A ■ B = ab cos 0 


(A7b) 


The cross product of two vectors, A x B , is a vector , C , whose 
components are defined by , 


A X B = C = 




- 1 
-Vs ) 

- a b 

2 1 J 


(A8) 
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It is plain that A x B = -B x A . 

The cross product of a vector into itself is the null vector 


A X A 



(A9) 


The operation of cross product and dot product commute, that is, 
given three vectors A , B and C 


(AxB)-C = A. (BxC) = B - (CxA) (AlO) 

It follows from ( AlO) that if C = A , or C = B , we have 


(AxB)-B = A* (BxB) = 0 
(AxB)*A=B-(AxA)=0 


(All) 


It follows from (All) that the cross product vector, A x B , is a 
vector perpendicular to both A and B. Thus A x B is a vector normal 
to the plane defined by A and B . 

The triple cross product of three vectors, (AxB)xC is a vector 
which is contained in the plane defined by A and B and is orthogonal to C. 
Thus, we have 


( A X B ) X C = B (A • C) - A (B • C) 


(A12) 


The square of the length of the cross product vector, A x B , is 
given by 


( A X B ) • ( A X B ) = (A • A) (B • B) ” (A • B) 


2 


(A13) 



I mil II 


Thus, we have 

I A X B I = a b I sin 0 I (A13a) 

where 0 is the angle between vectors A and B . 


B. Rotation of a Vector About a Un i t Ve cto r Through an Angle 

The differential equation for the instantaneous rotation of vector , R , 

A 

about a fixed unit vector , A , is given by 


“ R = w A X R 
dt 


(A14) 


where w is a scalar equal to the instantaneous rotation rate. 


To acquire the solution of Eq. ( A 14) we may interpret the operation 

A 

A X as a matrix. We have, from Eq. (AS), for any vector, B , 

“ -^3 

*3 ® '\) { ’’2 ) ( A15 ) 

A 

where a^ , a^ , and a^ are the components of A. 


A X B = 



A 


Thus, the operation , A x 

, is a 

skew S3Tnmetric 3x3 matrix given by 


/ 


\ 


1 0 

~a^ 


A 


3 

\ 

Ax = 

^3 

0 

-a^ (A15a) 


i -a„ 

a. 

0 


\ ^ 

1 

/ 

/ 


The matrix satisfies the characteristic equation 


A o 

( Ax)^ 


2 


- a 


(Ax) 


(A16) 
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To prove this, we note that 


A ^ 

Ax A X B 

Since 

A A 

A * ( A X B ) = 0 
and B is any vector, we have, for (Ax) the identity 



D a A A A A A 

= A A* (AxB)-A*A(AxB) 


(Ax)(Ax)(Ax) = - A-A(Ax) 


which proves Eq. (A 16). 


2 . 

Since A is of unit length, a =1, the eigenvalues of (Ax) 
0 , i ,-i ; where i = /^. 

The formal solution of Eq. (1. 14) is 


R(t) = 


w t ( A X ) 
e 


R<V 


Since the roots of the characteristic equation of (Ax) are 
i , and -i , we have, from the Lagrange interpolation formula, 


wt (Ax) _ 


= E 

k=l 


3 

n 

i = 1 [ (Ax) - Xjl ] 


3 

n 

i = i 


(X, -X.) 

k j' 


wt 


where 


I = 


(A17) 

(A17a) 

(AlTb) 

are 

(A18) 

X = 0 , 

(A19a) 

(A19b) 
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Evaluating Eq. (A 19) for X = 0 , i and -i , we have 

A 

^wt(A X) ^ ^ (A X ) (A X ) + sin wt A x ) 

A 

The inverse (or reverse) transformation is e since 

A A 

^-wt(A X) ^wt(A X) ^ ^o ^ j 

Furthermore, the inverse is the transpose since (Ax) is skew symmetric. 

It follows that the transformation is orthogonal and represents a rigid rotation 
of every vector R (t ) into a corresponding vector R (t) about the unit vector, 
A , through the angle, w t . 

If we choose the three unit vectors along the coordinate axes of the 

rotated system for R (t ) , the corresponding columns in the matrix, 
wt(A x) ^ 

e ^ , given by Eq. ( A 20) will be the images of the three axes in the fixed 

coordinate system. 

A 

Since the transformation depends only on wt and A we define 

A 

T ( a , A ) = (A22) 

A A 

If we have a sequence of rotations of cl^ about A^ , about A^ 
in the rotated system resulting from the first rotation, and ot about A 

A O 

in the system rotated about A , etc., we have 

A A A 

G = T ( , A^ ) T ( 0^ , A^ ) T ( , Ag) (A23) 

a transformation from the rotated system into the fixed original system. 

As an illustration of the method we choose an inertial fixed reference 
frame, given by 

x^ along the Earth North Pole; 

X at the intersection of the Greenwich Meridian and the equator at 
the initial time, t^ ; 

Xg at the intersection of the -90^ meridian and the equator at to . 


(A20) 


(A21) 
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* 


As the Earth rotates about its North Pole with constant rate, w 

we seek the transformation between the initial inertial reference frame and 

a rotating frame fixed in the rotating planet. We define the transformation 

to be , from Earth fixed to inertial. Let At = t - t 
IE o 

A 

GjE = At , x^ ) (A24) 


Evaluating Eq. (A20) with 



G 


IE 


0 cos w At 

£j 

1^0 sin Wg At 


0 

-sin At 
E 

cos w„ At 
E 


(A25) 


(A26a) 


It is convenient to define a level earth coordinate system with its 
X axis to the point on the Earth^s surface with right ascension, » 
and declination 6 . We obtain the transformation from level Earth to 
inertial coordinates through 


G 


IL 


= T (Qf^ x^ ) T ( 6^ ) 


(A26) 


In matrix form, Eq. (A26) is given by 


G 


IL 




cos 6 

sin 6 sina-j^ 
-sin 6 cos Qi ^ 


0 

cos a 
sin a 


sin 6 

-cos 6 sinog^ 
cos 6 cosc^l^ 


(A26a) 


In Earth level coordinates, the x„ axis is the radius vector from the 

3 

A 

Earth center in the inertial system; the x^ vector points due west and the 
x^ vector is due north. 
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C. Rotation Rates 


The matrix , T ( a , A ) , is an orthogonal transformation. Thus it 
leaves the length of every vector invariant. 


C R } = T { } 


(A27) 


then 

R • R = [ Rq ]^t"^T { = R^ • R^ (A27a) 

Furthermore, given any two vectors A and B , the angle between them 
0 , is invariant under the transformation of ( T ) . Let 

{ } = T { A 3 

and (A 28) 

{ } = T { B } 

then 

A^ * = a^ cos a (A28a) 

where a is the angle between A^ and B^. 

From Eq. (A 28) we have 

A^ • B^ = [ A ] t"^ T{B} = A- B = abcos0 (A28b) 

Since a^ = a and b^ = a we conclude that a = 0 . 

Since the cross product of two vectors is a vector, we have that 
(T) [AxB} = (T)[A}x(T)[B3 (A29) 
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Furthermore, we have 


(T){AxB} = (T)(Ax)t"^T{b}=(T){A}x(T)[B 


Since T [ B } is an arbitrary vector, we conclude that for all 
( T ) and { A } 


( T ) ( A X ) s ((T) C a } X ) 


The result enables us to relate the rotation rate in one coordinate 
system to the rotation rates in the transformed system. 

Let 


G 


T ( w^t 


Af ) T ( W2 t , A^ ) T ( Wg t , Ag ) 


We have that 
d G 

— ^ = w ( B X ) G 
d t 

A 

where B is the unit instantaneous rotation axis of the transformation, C 
and w is its scalar rotation rate. 

Differentiating Eq. (A31) we have 

J ^ A /V, A A 

G = (AjX) G + T (w^t.Aj) Wg (AgX) T (Wgt, A^) T (Wg t,A 

A A A A 

+ T (w^t,A^) T (Wgt.Ag) Wg (AgX) T (Wgt,Ag) 


} (A29a) 


(A30) 


(A31) 


(A32) 


,) 

(A33) 
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Equating Eq’s. (A32) and (A33) we have 


A A A A 1 ^ 

w (Bx) = w^(A^x) +T (w^t,A^) WgCAgX) t” (w^t,A^) 


(A33a) 


+ T(w^t,A^) T (W 2 t,A 2 ) T ^(Wgt^A^) T ^(w^t,A^) 


Using the vector cross product identity of Eq. (A30) we have 


A A ^ A A ^ 

w (B X) = w^(A^x) (w^t,A^) } X J 

+ Wg (^ { T (w^t, A^) T (Wgt, Ag) Ag } X ^ 


(A34) 


Since the identity is true for all vectors we may remove the cross 
product operator and obtain the desired instantaneous rotation rate vector 

A 

w B in terms of the components in the rotating systems. 


w B = 




T (w^t,A^) A^ + Wg T (w^t,A^) T (w^t, A^) A^ 


(A35) 


To obtain the inverse relationship we note that 

A = - w (B X) g“^ (A36) 

dt 


Following the same logic as in Eq’s. (A31) through (A35) we obtain 


w B = 


WgAg + Wg T (W2t, Ag) A^ + T (Wgt, A^) T (w^t.A^) A^ (A37) 


Equation (1. 37) may be used to obtain the body rotation rate in the 
local level coordinate system, from the knowledge of the measured 

body rotation rate with respect to the inertial system, . We have 


dt 


BI 




X T 


BI 


(A 38) 
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where 


^ A ^ AAA 

Tbi = T(i, -(p) T (), -0) T(k, - It) T(l, tr) Tfl, - 6) T(i . - a) (A39) 


It follows that 


Ogj = pi + 0 



r -sin 0 ^ 

( sin <;^cos 9 ) 
1 cos <p cos y 




+ a 



(A40) 


Since 


6 = sin 


-1 


and 


a = tan 


-1 


~y ^ 

zJ 


(A40a) 


We have 


6 = 


2 • 

X R-R - r X 

~2 2 2 i" 

r (y +z ) 


(A40b) 


a = 


- y - tan a z 

2 2 i 

(y +z ) 


We may now solve for the Euler angle rates, <p , $ , and 0, from the 
differential equation 


.A • r ^ ^ ' r ® 

<p i + 0 y cos 0 } + 0 ( sin ip cos 0 


'v. 


sin(p 


l^cos <p cos 0 




BG 



(A41) 
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APPENDIX B 


Waypoint Path Construction 

This appendix derives the equations for constructing a desired path 
consisting of straight lines (great circles on the surface of a spherical Earth) 
connected by arcs of circles, fixed on the surface of a rotating spherical 
Earth, and rotatii^ with the Earth. 

The input data required for generating the path is minimal, con- 
sisting of the latitude and. longitude of a sequence of way points, the radius of 
the connecting circular arcs, and the desired altitude and airspeed at each 
way point. 

The path is constructed by sequentially connecting each pair of way 
points by a great circle arc on the surface of the Earth. At each interior 
way point, a circle is constructed by the intersection of a circular cone from 
the center of the Earth to the surface of the Earth. The radius of the circular 
conical base is taken to be the input turn radius at that interior way point. The 
center of the cone is chosen so that the circle is tangent to both the incoming 
great circle and the outgoing great circle. 

The coordinates of the unit vectors from theEarth center to the two 
tangent points on the surface of the spherical Earth are designated as additional 
way points, and the unit vector to the middle of each turn replaces each interior 
way point. The vertical construction of the path is now accomplished by 
generating a piecewise linear altitude variation from the middle of each turn 
to the middle of the next turn as a linear function of the ground track distance 
between the middle turn points. The initial and final segments are the two 
great circles connecting the initial way point to the first incoming tangent 
way point, and the great circle connecting the outgoing final tangent way point 
to the final way point. 
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A. Initial Data 


Let N be the number of way points (N must be at least three). 
At each way point we must supply the following; 


X (I) = longitude in degrees 
6 (I) = latitude in degrees 
h (I) = altitude in meters 
Vj^I) = airspeed in meters/sec 

I == 1, N 

and at each interior way point we must supply the following; 


H^(I) 

I 


radius of the turn (meters) 
1, N - 2 


(Bl) 


(B2) 


B. V ector Representation of Each Way Point 


The vector representation of each way point is in the Earth fixed 
frame. Converting the angle data to radians, the way point unit vectors 
are given by, 


WR (I) = 


( sin 6 (I) 
^ -cos 6 (I) 
cos 6 (I) 


sin X 
cos X 



I = 1 , N 


(B3) 
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The Unit Normal V ectors 

Between each consecutive pair of way points we pass a plane which 
intersects the Earth surface in a great circle. The unit normal to that 
plane is given by 


WN (I) = 


WR (I) X WR (I + 1) 

■ ^ V. ' ‘ ■ 

|WR (I) X WR (I + 1) I 


I = 1, N - 1 


(B4) 


The vector magnitude of the cross product is required for the computation of 
the change in yaw heading angle at each interior corner way point. Let 


DN (I) = I WR (I) X WR (I + 1) I 
I = 1, N - 1 


(B4a) 


D. The Change in Yaw Heading 

The change in yaw heading at each interior way point is the angle through 

A 

which the unit normal vector WN (I) must be rotated about the unit way point, 

A A 

WR (I + 1) , in order to coincide with the outgoing unit normal vector, WN (I + 1), 

A 

(See Fig. 23). Since the interior waypoint vector, WR (I + 1) is peipendicular 
to both the incoming and the outgoing normals, we have 

A A A A 

WN (I + 1) = cos [ All) (I) ] WN(I) - sin [ Alb (I) ] WR(I+1) x WN(I) (B5) 
From Eq. ( B5) we have 

A A 

WN (I) • WN (I + 1) = cos [ Alb (I) ] (B6a) 

I = 1, N - 2 

and 

A A A 

WR (I + 1) X WN (I) • WN (I + 1) = ~ sin [A0 (I) ] (B6b) 
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Since, from Eq. (B4) 


WN (I) 


WR (I) X WR (I + 1) 
DN (I) 


we have 


(B6c) 


WR (I + 1) X WN (I) 


WR (I) - WR (I + 1)(WR (I) « WR (I + 1)) 
DN (I) 


and Eq. ( B6b) becomes 


WR (I) • WN (I + 1) 
DN (I) 


- sin I Ail) (I) ] 


1 = 1, N - 2 


(B6e) 


It follows that the change in heading is given by 



- WR (I) • WN (I + 1) 




1 = 1, N - 2 


(B7) 


The four quadrant definition of arc tangent should be used to evaluate Eq. (B 7). 


E. Center of the Circular Turn 


We note that for each interior way point WR (I + 1) , the unit radius 

A 

vector at the center of the turn CR (I) , the unit vector at the middle of the 

A 

turn WC (I) , and WR (I + 1) all lie in a plane. (See Fig. 24). Furthermore, 


A A A 

since WN (I) and WN (I + 1) are perpendicular to WR (I + 1), it follows that the 

A A A 


sum WN (I) + WN (I + 1) is perpendicular to WR (I + 1) and also lies in the 
same plane. 
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WC (I) 


Note: 



WN (I) + WN (I + 1) 


For a negative ( A0 ) the vectors WN (I) and WN (I + 1) will be in the 
opposite direction. 


Figure 24,— Center of Turn and Middle of Turn Vectors 
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We have 


A 

SUMV (I) 


A 


A 


WN (I) + WN (I + 1) 

A A 

|WN (I) + WN (I + 1) 


I = 1 , N - 2 


where 


A A 

WN (I) + WN (I + 1) 



(B8) 


(B8a) 


It follows that WN (I + 1) and SUMV (I) form a perpendicular basis for all 
vectors in the common plane. Let ^ be the angle between the center of the 

A A 

turn, CR (I) and WR (I + 1) ; then the unit vector at the center of the circle 
is 

A A A 

CR (I) = cos j8 WR (I + 1) - sign ( ^ ) sin jSSUMV (I) (B9) 


I = 1, N - 2 


Let be the angle between CR (I) and WC (I) on the circle with radius, 

Rj, (I). Then, we have 

Rp (I) 

sin = (BIO) 


and the unit vector to the middle of the turn, WC (I) , is given by 


WC (I) = cos ( jS - a ) WR (I + 1) 2- sign ( Aj/) ) sin ( jS - a ) SUMV (I) 

(Bll) 


1=1, N - 2 


To obtain an expression for the ar^le , /S , we have from spherical triai^les, 
using Fig. 23 , 


sin B 


sin ft 


cos 


M. 

2 


(B12a) 
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(B12b) 


cos j3 



i8 


From Eq. ( B12a) we conclude that since the maximum value of sin j3 = 1, , 
that the limiting value of is given by 


2 


-1 ^ 'N 

cos [sm ( — ) J 
e 


(B13) 


This value occurs when the two great circle paths are parallel meridians 
meeting at a point halfway around the Earth sphere from the circle center. 
For this limiting value of » the two radii from the center to the in- 

A /V 

coming and outgoing tangent points lie a great circle, containing CR (I) . 

The normal at the middle of the turn is required in the guidance 
parameter equations in order to compute the distance to go to complete the 
half turn. We have 


A A 

WNE (I) = SUMV (I) 
I = 1 , N - 2 


(B13a) 


For the end of the turn we note that the normal at the outgoing tangent way point 

A 


is WN (I + 1) . 


F. The Incoming and Outgoing Tangent Unit Vectors 

To obtain the equation for the unit vector tangent way points, we note 

A A 

that the circle center CR (I) and the normal to the great circle, WN (I) lie 

A 

in a plane containing the incoming tangent waypoint PI (I). (See Fig. 25). 

A A 

Since PI (I) is perpendicular to WN (I) , we have 


PI (21 - 1) 


1 r 

cos ^ tt WN 



I = 1 , N - 2 


(B14) 
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WN (I]^ 


Note: 



For negative ( A0 ) the directions of WN (I + 1) and WN (I) will be 
reverses in the figure. 


Figure 25. -Incoming and Outgoing Tangent Vectors 
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Similarly, for the outgoing unit tangent way point we have 


1 ^ ^ "N 

PI (21) = (CR (I) + sign ( Aj^)) sin aWN (I + 1) J (B15) 

I = 1, N - 2 

G. Altitude and Airspeed Gradients 

To compute the gradient in altitude and airspeed, we assign to the 
initial and final unit way points, as well as to each interior unit middle of 

A 

the turn way points, WC (I) , a desired altitude and airspeed. By deter- 
mining the ground track distance between these successive way points, we 
may compute the gradient in altitude and airspeed over each element. (See 
Figs. (26a) , (26b) and (26c). 

The distance between the initial way point and the first middle way 
point is given by 

DST (1) = r^sln"^ (| WB (1) x PI (1) |) + (1) I (B16) 

The distance between each of two successive middle of the turn unit way 
points is given by 

DST (I) = R^(I-1> |-|^(I-l)| + rgSln"’^(|pi(2I-2)xP(2I-l)|) 

+ (I) I (I) 1 (B17a) 

I = 2, N - 2 
(N> 3 ) 

Finally, the distance between the last middle of the turn way point 

A 

and the final touch down way point WR (N) , is given by 

DST (N - 1) = (N - 2) I -^^(N - 2) | + r^ sin"^ Q PI (2N-4) x WR (N) Q 

(B17b) 
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The gradient in altitude for each element is given by 


grad (I) = . ML l l l -hj l) 
DST (I) 

I = 1, N - 1 


(B18) 


The gradient in airspeed is given by 


DVG(I) = 


Vd(I^1)-Vd(I) 

DST (I) 


(B19) 


H. Way Point Guidance Array 

If the initial input contains N way points, the final path contains 3N - 4 
way points. These way points subtend 3N - 5 segments; N - 1 of these are 
arcs of great circles and 2 (N - 2) comprise N - 2 pairs of half turns. Each 
segment is described by a vector array of 19 elements used in RNAV 
guidance on that particular segment. 

The guidance parameters, detailed in Section I (b) (page 12) of this report, 
require data which are constant over each of the 3N - 5 segments. These data 
are conveniently arranged in 3N - 5 guidance parameter vector arrays to be 
computed initially at the beginning of each flight and called by the guidance 
logic in sequence as required. The elements of each array are listed below 
for both great circles and half turn circles. 

Let 1 1 go from 1 to 3N - 5 

J = Integer part ( 1 1/3) + 1 
K = 2J - 1 

If MODULO (I 1, 3) = 2, then L = J ^ 2 Q) 

If MODULO (11,3) = 0, then L = J - 1 
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For Great Circles 


WP (1) = 

WN (J, 1) 

WP (2) = 

WN (J, 2) 

A 

WP (3) = 

WN (J, 3) 

A 

WP (4) = 

PI (K, 1) 

A 

WP (5) = 

PI (K, 2) 

A 

WP (6) = 

PI (K, 3) 

= 3N - 5, 

(4) = 

WR (N, 1) 

A 

WP (5) = 

WR (N, 2) 

A 

\WP (6) = 

WR (N, 3) 
y 

WP (7) = 

sign (Aj/) 1 

WP (8) = 

CR (J, 1) 

A 

WP (9) = 

CR (J, 2) 


(B20) 

Cont’d. 


WP (10) = CR (J, 3) 


WP (11) 
WP (12) = 
WP (13) = 


- 0 


(J) 

0 


-1 


WP (14) = sin (tan (GRAD(J) ) 



H (1) 

VEND = (1) 

D ^ 


-1 


If 1 1 > 1 , but < 3N - 5 

/^END = HEND + GRAD(J)(r sin*'^ (| PI(2J-2) x PI(2J-1)|) 
\^END = VEND + DVD(J) sin' (| PI(2J-2) x PI(2J-1 )|) ) 


If I I = 3N - 5 

/^END = H (N) 
V^END = V^ (N) 



WP (15) 

= HEND 

WP (16) 

= GRAD(J) 

WP (17) 

= VEND 

WP (18) 

= DVD(J) 

WP (19) 

= J 

WP (20) 

= |Ai1)(J>/2| 


For Half Turn Circles 
If MODULO ( 1 1,3 ) = 2, then L = J 


WP (1) 

WP (2) 

WP (3) 

If MODULO (I 1,3 ) 

WP (1) 

WP (2) 

WP (3) 

WP (4) 

WP (5) 

WP (6) 

WP (7) 

WP (8) 

WP (9) 

WP (10) = 

WP (11) = 
WP (12) = 

WP (13) = 

WP (14) = 

HEND 
VEND 


WNE (L, 1) 

A 

WNE (L, 2) 

A 

WNE (L, 3) 

= 0, then L = J-1 

A 

WN (J, 1) 

A 

WN (J, 2) 

A 

WN (J, 3) 


0 

0 

0 

sign (Aij) (L) ) 

A 

CR (L, 1) 

A 

CR (L, 2) 

A 

CR (L, 3) 

I Al/) (L) / 2 I 


\ (L) 

0 

sin ^an"^ (GRAD (J ) )) 

HEND + GRAD (J) ^^(L) | A0 (L) / 2 | 
VEND + DVD (J) ^R,p(L)| A0 (L) / 2 


(B20) 

Cont’d. 


(B21) 
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I mmi ■■■■Hill 


WP (15) = HEND 

WP (16) = GRAD (J) 

WP (17) = VEND 

WP (18) = DVD (J) 

WP (19) = J 

WP (20) = 0 


(B21) 

Cont’d. 
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APPENDIX C 


Discrete Formulation of a Third 
Order Complementary Filter 

This appendix derives a closed form solution of the third order comple- 
mentary filter for use in implementing a discrete formulation of the filter for 
flight computers. The solution permits the use of input gains to yield either 
3 real or 1 real and 2 complex conjugate roots. It provides for fixed coef- 
ficients to update the aircraft state position, velocity and acceleration directly 
from input noisy accelerometer data, and navigation aids such as MLS, VORTAC, 
and barometric data, etc. , without requiring Runge Kutta integration. 

In what follows, the measurements are converted into a runway Car- 
tesian vector ( X, y, z ). The details of the transformation are contained in 
Section II of this report and is not repeated here. Each component of the 
coordinate system is estimated by means of three first order linear differen- 
tial equations driven by the converted measurement, say x , and the measure- 
ment of acceleration in that coordinate, *x* . 

A . The Thi ^ O rder Complem enta ry Filter 

The differential equations for the third order complementary filter for 
one of the three coordinates, (say x (t) ) may be written as 


^2 k :2 ^ ^ ^ ^ ^3 

X3 = Kg ( X (t) - x^ ) 

where 

Kf , Kg and Kg are the filter gains. 


(Cla) 

(Clb) 

(Clc) 
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X (t) is the navigation aid pseudo measurement of the coordinate, in 
the X direction. 

X (t) is the measured noisy acceleration, assumed constant over the 
internal 0 ^ r . 

From Eq. ( Clc) we may solve for x^ - x (t) 


- X <t) == - 



(C2) 


Differentiating Eq. (C2) we get 


Xi-x(t) = - K- Xg 

o 


(C2a) 


Substituting for from Eq. ( C la) we obtain 


(X (t) - x^ ) + x^ - X (t) = - — Xg 


(C2b) 


We may replace (x (t) - x^ ) in Eq. ( C2b) by its equivalent in Eq. 
( C 2) and obtain 


’'2'’=^ = - IT 

o 


(C3) 


Differentiating Eq. (C 3) we obtain 


X2-x(t) = - — (x-^K^Xg) 
o 


(C3a) 


Replacing x by Eq. (C lb) , we obtain 


Kg ( X (t) - x^) + Xg + X (t) - X = - — ( Xg + Xg ) (C3b) 

3 

Again we replace x (t) - x^ in Eq. (C 3b) by its equivalent in Eq. (C 2) and 
obtain 


X(t)-X(t) = (Xg+K^ ^^3 ^ Kg ^3 ^ *^3 ^3 ) 

o 


(C4) 
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We now make the assumption that the navigation aid measurement 
X (t) is modeled as a quadratic polynomial of time. 

. . . .2 

x(t)=x + x t+x (C5) 

o o o 2 

We determine the initial coefficient as follows: 

At t = 0 (the beginning of the discrete time interval). 

X (0) = = Xg (0) 

• • • • 

X (0) = x^ = X ( r ) + Xg (0) 

At t = T we set the navigation aid measurement equal to x ( r ) • 

2 

X ( T ) = x^ + Xg (0) T + Qx ( T ) + x^ (0) J (C5c) 

We may now solve for the initial value of x (0) , namely x^ 

2 

x^ = X ( T ) - (0) T - (y ( T ) + Xg (0) ^ (C6) 

Subtracting x^ (0) from both sides we have 

2 

Xo - Xf (0) = X ( T ) - (xj^ (0) + Xg (0) r + ( Xg (0) + X ) ^ (C6a) 

Thus we have succeeded to model the observation x ( r ) in such a 
manner that if the difference between the observation x ( t ) and its pre- 
diction, as a quadratic function of x^ (0) , x^ (0) , x^ (0) and x ( t ) 
vanishes, then the initial value of x^ - x^ (0) also vanishes. 


(C5a) 

(C5b) 
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Returning to Eq. (C4), we note that the right hand side is a constant 
and using Eq. (C 5b) , Eq. (C4) becomes 


■^3 ^ = S^3<«> 


(C4a) 


The solution of Eq. (C4a) is either 


X - X (0) = a e ^ ^ + b e ^ cos 5t + c e ^2 ^ ^ (C7) 

o o p 


(if the gains K , K , and K produce 1 real root and a pair of complex 

X ^ O 

conjugate roots), or else 


X - X (0) = a e + b e + c e 

O O 


(C8) 


if K- , K„ and K„ produce 3 real roots of the characteristic equation 

d. ^ o 


+ K + K z + K = 0 

X ^ O 


(C9) 


For the case of only one real root, Eq. ( C7), we have 


( z )(^( z +Qt^ )^ = 0 


(C7a) 


and 


Ki = «i + 2 <*2 


"^2 = ^“ 2*1 ® 2 ^ 


Kg = 


(C7b) 

(C7c) 

(C7d) 
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For the case of 3 real roots, Eq, (C 8 ) , we have 


(z+Og) ( z + Ofg ) = 0 (C 8 a) 

Oil + «2 Kf (C 8 b) 

«1 «2 ^ ^ ®2®^3 ~ ^2 (C 8 c) 


We now proceed to solve the discrete third order complementary filter 
for the two cases. 


Let the roots of the characteristic equation be a. , cv„ + i ]8 , 
Odo - i (3 ; then, the solution of Eq. (C 4a) is Eq. (C 7). We may now 
differentiate Eq. (C7) to obtain a solution of Eq. (C2), We have 


K 3 ( Xf - X (t) = X 3 = + b e ^ cos /3 t + sin |3 t ) 


a e 


«2 


a. t 
1 




(CIO) 






, _ , sin B t _ , ^ 

+ c e ( a 2 — cos jSt ) 


Similarly, we may obtain an e 3 q>ression for x + K x . Recalling that 

O X O 

K, = a. + 2 ao we may solve Eq. (C 3) and obtain 
J. X ^ 


r -°i* 

be ^ 


K3(x2-5(t))=-X3-K^X3=4b e ^ a^O^) oos ^ 

+'c e -a^cos/3t 


>(C11) 
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We may combine Eq’s. (CIO) , ( Cll) and (C7) in matrix form to 
obtain 


( X ^ - X ( t ) P 


r 


K (x^ - X (t) ) 


a e 


-«it ■>! 


}= (A) <b 


” « t 

2 


^.=^3 - =^3 y 


c e 


- t 


where the elements of the 3x3 matrix (A) are given by 


‘11 


‘12 


13 


21 


22 


23 


31 


32 


33 




Oi cos t + 3 sin ^ t 


sin 8t 

a 2 — — - cos 8 t 


2 a , a « 

12 


2 2 

i & 2 +a^a2)-cosj3t + o' 


,22 ^ sin R t 

< «2 ^ “l ®2> "It " ®1 


^ i3 sin j3 t 
cos jS t 


cos jS t 
sin jg t 


(C12) 


(C12a) 
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Setting t = 0 , Eq. (Cl2) becomes 


j 

«1 

«2 

-1 1 

i a 

1 

1 

0 i 

! 


0^2^ + 

-“i 

1 

<'■ b 
1 1 

L« J 

1 

1 

"1 

o 

o 




(C13) 


The solution for a , b, and c is given by 
* = K3 (x^(0)-5^) 

b = K ( X (0) - X ) (C14) 

2 2 
-ot^ - B 

c K 3 (x^< 0 )-x^) 

where A is the determinant of the matrix (A) 

A = - 2 a 2 (C14a) 

From Eq. (C 6a) we may substitute for x^(0) - in Eq. ( C14) and obtain 
a = Kg ( x^ ( T ) - X ( T ) ) 

b = -^K3 (x^(t)-x(t)) <C14b) 

2 2 
cuot.-Oio-B 

c = Kg ( x^ ( T ) - X ( T ) ) 

2 

Xf = x^ (0) +Xg (0) T + ( Xg (0) +x ( T ) ) — (C14c) 
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After substituting the solution for a , b , and c into Eq. ( Cl2) we 


{ (0) + (Xg (0) + X ) T > + ^ (x^ ( T ) - X ( T ) 




X3(0) 


Gg (Xi ( T ) - X ( T ) 


where 


/ 2 2 ~^2 ^ \ 
I e -20^0^2)^ COS/ 3 T ' 




Z' 2 2 2 2 ^ 2 sin/?T ' 

+ («j (Og -S )- (0(2 )a2) e — j- .- 


/ 2 “°i "®^2 

I 2 Og ( e - e COS 5r ) 


\ /^ 2 2 2 2 : 
\+Qo(i (0(2 -f)-(0l2 +$ 


, 2 , ^ “^^'’’sin^T 

) Z® -r- 


-a^T 

O'- ( e - e cos )3 r ) 


1 2 2 

\ +ag “ ^ ® 


si 


sin gT A 
B 


A somewhat more convenient form for the solution is obtained by replacing 
by G^ - 1 . Then Eq. ( C15) becomes 


(T ) 


< X2(T) 


"I f X ( T ) 

1 ^ 

1 


(G^-1) (x^ (r) - X (T)) 

1 

] 

> = ^ ^2 (0) + (XsCO) +x ) T 


S ^ ^ 

i 

J Us <*» 


A 

^Gg (Xf (T) -X (T)) ^ 
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In this form it is clear that if the residual (t) ~ x (t) is zero, 
then the solution of (t), ( r) and x^ (t) is the so called "dead 

reckoning" solution obtained for a linear system under constant acceleration. 
Let the roots of the characteristic equation be a- , and a . 

J. ^ O 

In place of Eq. ( CIO) we now have 

^ 

Kg(xi - X ( r)J = a e 

In place of Eq. ( Cll) we now have 


. "0^1 ““2*^ 
Kg(x2 (r)} = (02 +ag )'e +bo2 («i +« 3 ) ® 

-Q?3 T 

+ c q?3 ( 0^^ +®2 ) ® (C18) 

For the equation for x - x (0) we have from Eq. ( C8) 

o 

-O^T ~^2'^ 

X - X (0) = a e +b e + c e (C19) 

o O 

The matrix equation replacing Eq. ( C12) is given by 



- X (T) y 


r 

®2 

A 

^3 


r 

a e 

Is 

(x2-J (T)) 


Oi (Og +«3> 


«3 <“l ■'“2 

V 

b e 


- ^3 (“> 



1 

1 

1 ~QL T 

3 ^ 

!c e 

L 


(C20) 


Setting t = 


a = 


0 in Eq. ( C20) and solving for a , b and c we obtain 
Oi 

K^(x,(0)-i ) 




3 ' 1 


Of.. 




(C21) 


Ur 


c = 


K (X (0) - X ) 

(«3-ai) (Qig-C^) o' 
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The final solution equivalent to Eq. (C 16) is given by the G. 


coefficients 
2 


Gl-1 = 


a. 


“ttlT 2 
e 1'^ 0 ^ 




-OfiT 2 

e 0^ 

+ 


“oe T' 


(ag-ai) (c^^2) 


-1 


e 02 ^ e“^2T a^(Ci^+a^) e “2'’’ ^ 

°2 " (a^-a^) (a^-Og) (ag-a^) (o^-o^) ^ 


G. = 






e ^ K 

+ 


.“0^3 T 




B. Determination of the Characteristic Roots From the Filter Gains 


It is sometimes more convenient to input the three gains, K. , rather 
than the roots of the characteristic equation. In that case it is necessary to 
determine the values of the three roots and to find whether two are complex. 
Let the gains be K , K and K , The characteristic equation is 

J . M O 


3 2 

z + K. z + z + =0 

12 3 


(C23) 


with 


K, s 0 
1 


(C23a) 


First form the discriminant to test for reality of the roots. 


Let 

C = (C24) 

P = *^ 2 - g /3 

Q = C) Ki 
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The discriminant is given by 


D 


3 2 

4 P + 27 Q 


If D < 0 3 real negative roots 

If D ^ 0 1 real negative and 2 complex conjugate 


For D < 0 



For D i 0 
B 
E 
F 
G 
H 



D 

^ 108 ^ 


A cos 9 

“ { cos 6 sin 6 

(cos 9- /"s’ sin 9 

2 


Q +B 


- i Q - B 


sgn (E) 


1 

3 


sgn (F) 



1 

T 


(C25) 


(C26a) 

(C26b) 


(C27) 


(C28a) 

(C28b) 

(C28c) 

<C28d) 

(C28e) 
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iri ri 


1 

'^2 = — ^ -^(G+H) 




2 




(C29) 
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APPENDIX D 


Wind Estimator Models 


The physical model used to estimate the winds is based on the assump- 
tion that the aircraft attitude is maintained in a trimmed condition, aligned 
with the relative wind. In this equilibrium condition, we have for the horizontal 
components of the wind 


w 

X 




V . J cos 0 
airspeed ^ 


(Dl) 


w 

y 


Rg(2) 


V . , sin 0 

airspeed 


(D2) 


If we retain the initial value of w and w as the best estimate of the con- 

X y 

stant winds in a local level coordinate frame, we may model the deviation from 
this original estimate as an exponentially correlated random variable. 


Thus for the Kalman filter we have 

bw (t) = bw ( t - At) e ^ +0'. 

X X 

At 

bw (t) = bw (t-At)e ^ + o / ^ 

y y t 



where 


(D3) 


(D4) 


cr =1. m/sec , t = 100. sec , At = .05 sec 
w w 


The best estimate of the horizontal winds for the Kalman filter is 


Wx (t) = (to) + bw^ (t) (D5) 


Wy (t) = Wy (to) +bWy (t) (D6) 
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The observation of the winds are supplied by Eq*s. ( Dl) and (D 2) and 
the residual for the Kalman filter equations are given by 


Aw 

X 


w (y) 

X 


Eq. (4. 1) 


W^<t) 


(D7) 


Aw 

y 


w (t) - w (t) 

y Eq. (4. 2) y 


(D8) 


The Kalman gains are computed in a manner already described in Ref. 1 (page 31) 
and are not repeated here. The updated value of the wind biases are given by 


A A 


bw (t) 
x' ' 

= bw (t) + A bw 

x' ’ X 

(D9) 

A 

bw (t) 

y' 

A 

= bw (t) + Ab w 

y y 

(DIO) 


For the case of the complementary filter we chose a first order ejq)on- 
ential model without reference to any estimate of the variance. We have 

- At/t^ 

b w (t) = bw (t - At) + e ^ (w (t^) + bw (t - At) - w (t) (DH-) 

^ ^ XX ^Eq. (Dl) 


- At/t„, 

bw (t) = bw (t - At) + e ^ (w (t^) + bw (t - At) - w (t) (D12) 

^ ^ y y yj.q_ 2) 


where bw^ (t^) = bw^ (t^) = 0 , 


tw = 100 sec, At 


. 05 sec. 
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